vdW-DF study of energetic, structural, and 
vibrational properties of small water clusters and ice Ih 
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We present results for a density functional theory study of small water clusters and hexagonal 
ice Ih, using the van der Waals density functional (vdW-DF). In particular, we examine energetic, 
structural, and vibrational properties of these systems. Our results exhibit excellent agreement 
with both experiment and quantum-chemistry calculations and show a systematic and consistent 
improvement over standard exchange-correlation functionals — making vdW-DF a promising candi- 
date for resolving longstanding difficulties with density functional theory in describing water. In 
addition, by comparing our vdW-DF results with quantum-chemistry calculations and other stan- 
dard exchange-correlation functionals, we shed light on the question of why standard functionals 
often fail to describe these systems accurately. 



I. INTRODUCTION 



Despite the overwhelming abundance of water on Earth's 
surface and its composing the majority of all known life, 
much of the physics of water remains a mystery. Un- 
raveling that mystery requires the ability to take mi- 
croscopic information and extend it toward the macro- 
scopic regime. Classical force fields and other parame- 
terized models can calculate average macroscopic prop- 
erties of water wellji - — but generally yield little insight 
into the microscopic physics underpinning these proper- 
ties. On the other hand, quantum-chemistry methods 
such as M0ller-Plesset perturbation theory and coupled- 
cluster techniques have been invaluable in understanding 
numerous molecular systems and indeed are responsible 
for much of what is understood about water todayj^Z 
But, their poor scaling with system size limits their use to 
systems containing only a few molecules. 8 Density func- 
tional theory (DFT) scales substantially better and is 
thus, in principle, well suited to make the connection to 
the macroscopic regime. Unfortunately, standard func- 
tionals fail to adequately describe the weak van der Waals 
interactions that are critical to obtain accurate diffusive 
and structural properties of waters- 
Many of water's interesting properties stem from its 
ability to form complex hydrogen-bonded structures, the 
description of which requires an accurate treatment of 
dispersion interactions. Historically, DFT has done a 
poor job of treating these weak interactions. Recent in- 
terest in van der Waals systems has led to the develop- 
ment of a number of successful approaches aimed at im- 
proving DFT's ability to incorporate dispersion interac- 
tions accurately. These include post-processin g 14 i 15 and 
semi-empirical method s 16 ! 17 as well as hybrid 18,19 and 
non-local functiona l 20 ! 21 Here, we focus on the recently 
developed van der Waals density functional (vdW-DF), a 
truly non-local functional which fits seamlessly into DFT. 
This functional has been successfully applied to a variety 
of other weakly binding systems^ 2 . - — and is known to 
hold promise for improving DFT's description of the hy- 
drogen bonding in water *22r— Here, we use it to perform 



systematic calculations on small water clusters (H2 0) n 
with n = 1 — 5 and standard ice Ih . Comparison with ex- 
periment and second-order M0ller-Plesset perturbation 
theory (MP2) calculations demonstrates that DFT, uti- 
lizing vdW-DF, is able to obtain excellent results for a 
wide variety of water's properties, even in large, bulk-like 
systems. 



II. COMPUTATIONAL DETAILS 

Our DFT calculations were carried out using the plane- 
wave self-consistent field (PWscf) code within the 
Quantum-Espresso package^ utilizing ultrasoft pseu- 
dopotentials. For comparison, in addition to vdW-DF 
we used a standard local functional (LDA) and a semi- 
local functional (PBE). The three functionals all used 
Slater exchange and Perdew-Wang correlation 36 to de- 
scribe local exchange and correlation. For the PBE 
functional, the gradient correction approach of Perdew, 
Burke, and Ernzerhof was used. 37 For vdW-DF, semi- 
local exchange was provided by a revised PBE scheme 38 
and non-local correlation was provided by vdW-DF j 2Q i 21 
Wavefunction and charge-density cutoffs were 35 Ry and 
420 Ry, respectively. A self-consistency convergence cri- 
terion of at least 1 x 10 -10 Ry was used for all water 
cluster calculations and 1 x 10 -8 Ry for ice Ih- Ini- 
tial structures were relaxed until all force components 
were less than 1 x 10 -5 Ry/a.u. in the water clusters and 
1 x 10 -4 Ry/a.u. for ice. Stress relaxations were carried 
out in all ice calculations by relaxing each free lattice 
parameter until all stress components were smaller than 
1 kbar. To minimize interactions between periodic im- 
ages, simulation cells for the finite systems were sized to 
ensure a minimum separation of 10 A between atoms in 
one cell and atoms in neighboring cells. 

MP2 calculations were carried out using the quantum- 
chemistry package Gaussian^ employing Dunning's 
augmented, triple-zeta basis set (aug-cc-pVTZ). It has 
been shown that this basis set is sufficient to yield ac- 
curate properties for small water clusters, 4 and simple 
testing with the larger aug-cc-pVQZ basis set confirmed 
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only minimal changes in the properties of interest here. 
Additionally, the MP2 binding energies obtained here are 
within 5% of those reported by Santra et al4£ who ex- 
trapolated to the complete basis set. While counterpoise 
corrections were deemed negligible for most properties 
of interest, we found that they alter the vibrational fre- 
quencies by as much as 6% in the water dimer. For this 
reason, counterpoise corrections were applied to all MP2 
calculations carried out in this work. Initial water clus- 
ter structures were relaxed until all forces were less than 
4 x 10 -6 Ry/a.u. Detailed positions for all optimized 
structures are provided in the supplementary materials. 

Vibrational frequencies for the MP2 method were cal- 
culated directly, using Gaussian. For the DFT calcula- 
tions, forces on all atoms were calculated when each atom 
was displaced by ±0.0025, ±0.005, ±0.0075, and ±0.01 A 
along each of the three Cartesian directions. The dynam- 
ical matrix was then built by taking the first derivative 
of the forces using a nine-point numerical derivative. All 
systems stayed within the harmonic regime for all dis- 
placements employed. 



III. RESULTS 



Binding, Structure, and Electric Dipole 
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FIG. 1. (a) Binding energy per hydrogen bond for each of the 
multimers and ice Iu- The horizontal black lines correspond 
to MP2 calculations for the clusters and experiment for ice.— 
(b) Average oxygen-oxygen distance in each system. The 
black boxes represent experimental values . 42 ' 43 (c) Percent 
error relative to MP2 calculations of the dipole moment. The 
dipole moment of the tetramer is not shown since it is zero 
by the S4 symmetry of the complex. 



We start by presenting results for the binding energies of 
the water clusters (H2 0) n with n = 2 — 5, and standard, 
hexagonal ice Ih- The results are shown in Fig. QJa), 
along with MP2 values for the water clusters and the ex- 
perimental value for ice Ih- A1 Also shown for reference 
are values from two of the most popular local (LDA) and 
semi-local (PBE) functionals in use today. As can be 
seen in the figure, the LDA and PBE functionals exhibit 
a distinct overbinding. This is also evident in Fig. DJb), 
which shows the average oxygen-oxygen distance com- 
pared with experiment for both the clusters 42 and ice 
7^4^ The vdW-DF values are in excellent agreement with 
experiment and show a clear improvement over LDA and 
PBE, which draw the oxygen atoms too close together, 
consistent with the overbinding seen in Fig. QJa). 

Interestingly, a recent study by Wang et al.— has 
pointed out that, despite showing a substantial improve- 
ment over other functionals when calculating the self- 
diffusion and density of liquid water, vdW-DF tends to 
understructure the liquid. This understructuring mani- 
fests itself as a lowered density in the second coordination 
shell of the oxygen-oxygen radial distribution function. 
Wang et al. point out that this understructuring is likely 
an artifact of the chosen semi-local piece of the exchange 
functional, rather than a problem within the non-local 
correlation functional of vdW-DF. Many groups are in- 
vestigating the effects of different exchange functionals 
on vdW-DF and it will be interesting to see if the under- 
structuring problem can be solved. 

Another important quantity is the bulk modulus — a 



property that depends on the curvature of the energy 
surface when the system is subjected to isotropic expan- 
sion or contraction. In a recent study, 29 vdW-DF was 
used to calculate the bulk modulus of ice 1^. This study 
found vdW-DF to be in good agreement with the exper- 
imental value and serves as an independent extension of 
the findings presented here. 

To probe the accuracy of the electronic structure itself, 
the dipole moments were calculated for the finite systems. 
Figure QJc) shows the errors obtained (relative to MP2 
calculations) for the three functionals studied here. Ex- 
perimental results could not be found for clusters larger 
than n = 2, 44 but MP2 results agree with available exper- 
imental results to within 1%, so MP2 is used as the refer- 
ence throughout Fig. 0Jc). Detailed values for all dipole 
moments are presented in the supplementary materials. 
Again, results obtained from vdW-DF are in excellent 
agreement with these high-level quantum- chemistry cal- 
culations, while LDA and PBE show substantial errors. 

One can understand the better performance of vdW- 
DF since it is a non-local functional, but the physical 
ramifications of this distinction can be seen in Fig. [2] 
The figure shows the calculated charge density along a 
line from oxygen to hydrogen through both the cova- 
lent bond (solid lines) and the hydrogen bond (dashed 
lines). Not surprising, the charge densities corresponding 
to the covalent bond for each functional are nearly indis- 
tinguishable. However, when looking at the curves cor- 
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FIG. 2. Calculated valence charge density for the water dimer 
along the lines between the indicated oxygen and hydrogens 
through both the covalent bond (solid lines) and the hydrogen 
bond (dashed lines) for the three different functionals. 



responding to the charge density in the hydrogen bond, 
the cause for the overbinding of LDA and PBE — evident 
throughout Fig. [T| — is revealed. Both PBE, and to an 
even greater extent LDA, shift the peak of charge den- 
sity from the oxygen toward the hydrogen, resulting in 
an increase in the covalent character of the bond. (Note 
that, in the critical region between 20% and 60% along 
the path, LDA overestimates the charge density by as 
much as a factor of 2.) This, in turn, strengthens the 
bond and creates an overbinding that permeates all re- 
sults calculated with these functionals. 



B. Vibrational Frequencies 

Sit and Marzarii^ calculated the vibrational frequencies 
of the water dimer using the PBE functional. They found 
that PBE performed reasonably well, but significant dis- 
crepancies from experiment were found for some modes, 
particularly low- frequency, inter- molecular modes. They 
postulated, as have others, 32 that the problem arose at 
least in part from the inability of local and semi-local 
functionals to correctly treat the hydrogen bond. To ex- 
tend this investigation, we performed similar frequency 
calculations for all the water clusters — monomer through 
pentamer — using all three functionals and MP2. (Our re- 
sults for the dimer using PBE are very similar to those 
of Sit and Marzari.) The results for the larger wa- 
ter clusters exhibit trends similar to those found in the 
dimer. Figure [3] shows our DFT-calculated and avail- 
able experimentally-determined vibrational frequencies 
for the monomer^ dimer ^ and trimer^ reliable ex- 
perimental results for larger clusters could not be found. 
Results for all clusters compared with MP2 and the fre- 
quencies of ice Ih calculated with the three functionals 
are collected in the supplementary materials. It should 
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FIG. 3. Vibrational frequencies calculated with the LDA, 
PBE, and vdW-DF functionals compared to experimental val- 
ues for the water monomer,— dimer,— and t rimer.— The ver- 
tical line denotes a separation between inter- molecular domi- 
nated and intra- molecular dominated modes. The region be- 
tween cm -1 and 200 cm -1 (shaded area) was not included 
in the analysis. A detailed analysis of the dimer modes la- 
beled mode a, f3, 7, and S is given in Figs. 2] and [5] A full 
tabulation of mode frequencies for all methods can be found 
in the supplementary material. 



be noted that the frequencies calculated in this work were 
obtained using the harmonic approximation and no cor- 
rections were made for anharmonic effects. Such effects 
are known to be present in experimental results and can 
be influential in setting the precise frequency of various 
oscillations. Calculations including such effects have re- 
cently been reported at the PBE and PBE0 level^ but a 
direct comparison is not possible since the frequencies re- 
ported are for the deuterated monomer and dimer. Typ- 
ically, the changes induced by anharmonic effects are on 
the order of a few percent and fair comparisons can still 
be made to experiment. Nevertheless, in this work, truly 
quantitative comparisons of DFT results were done with 
respect to the harmonic approximation within MP2. 

As can be seen in Fig. [3j all DFT methods get simi- 
lar results for the monomer frequencies — not surprising 
since this is a system consisting solely of covalent bonds. 
Significant deviations occur in some modes for the larger 
clusters, however. These deviations are worse in the low- 
frequency region where inter-molecular interactions dom- 
inate. Most of the intra- molecular-dominated frequen- 
cies are obtained satisfactorily with all three functionals, 
but a few show significant error in LDA and PBE. Even 
within the intra-molecular modes there are oscillations 
which substantially change the hydrogen-bond geometry, 
so a proper treatment of non-local interactions is neces- 
sary to obtain quantitatively correct frequencies. Overall, 
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FIG. 4. Energy changes according to Eq. JTJ) as the water 
dimer is forced along the normal modes marked a and 7 in 
Fig. [3] The x-axes are on the same scale and represent a 
dimensionless progress variable along the mode. The insets 
above show the motion of the water molecules for each partic- 
ular mode. The grids lie in a plane containing the hydrogen 
bonding O-H-0 in the equilibrium position. In this figure, 
the hydrogen donor molecule is designated Ml and the hy- 
drogen acceptor molecule is designated M2. 



vdW-DF shows a consistent improvement for a majority 
of the modes. A full tabulation of modes for all methods 
and systems (including ice Ih) is given in the supplemen- 
tary materials. 

While some modes are well reproduced with all func- 
tional, others show a distinct spread in frequencies for 
different functionals. This is to be expected for inter- 
molecular modes, but it is somewhat surprising for the 
intra- molecular modes. This spread is caused largely 
by the different functionals' varying ability to accurately 
represent the hydrogen bond, as can be seen by a deeper 
analysis of some representative modes of the dimer sys- 
tem. Four such modes — labeled mode a, mode /?, mode 
7, and mode 8 — are marked in Fig. [3j Modes a and 7 
show a large spread in frequency for different function- 
als, while modes /? and 8 exhibit much less variation. We 



FIG. 5. Same as Fig. [4) except here shown for the modes 
marked /3 and 8 in Fig. [3] 

analyze these modes in greater detail as follows. 

The change in energy of any dimer configuration rela- 
tive to the equilibrium configuration can be written as 

AJStotai = A£ M 1 + AE M 2 + A£ int , (1) 

where A£mi and AEm2 are the strain energies within 
monomers 1 and 2, respectively, and AEjnt is the dif- 
ference in their interaction energy. As the atoms are 
moved along a certain mode, the individual contributions 
in Eq. (pQ) can increase, decrease, or remain unchanged, 
depending on whether they enhance, inhibit, or are ir- 
relevant to the change in energy for small excitations of 
that particular mode. The change in energy coming from 
these contributions relative to the overall energy change 
determines how much each energy term contributes to 
the vibrational frequency. 

Figures 3] and [5] show the contributions of the vari- 
ous energy terms in Eq. (pQ) as the dimer is displaced in 
both directions along the normal modes a, /3, 7, and 8. 
This was done in step sizes of 0.01 times the (normalized) 
eigenvectors, to a maximum value of 0.1 times the eigen- 
vectors. Mode a is comprised mainly of an out-of-plane 
oscillatory motion of the hydrogen donor water molecule. 
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Oscillations along this mode clearly change the geometry 
of the hydrogen bond. Since the overall energy change 
upon movement along this mode is relatively small, the 
interaction energy makes up a significant portion of the 
total restoring force. As the different functionals have 
a varying degree of success in describing the interaction 
energy, a significant spread in E'totai an d the frequencies 
results for different functionals. The motion for mode 
7 consists mainly of a symmetric stretch of the hydro- 
gen donor molecule. Again, this substantially changes 
the interaction energy as the dimer oscillates along the 
mode, resulting in a spread in Etotai an d the frequen- 
cies. On the other hand, mode j3 in Fig. [5] is comprised 
mostly of an angle flex in the hydrogen donor molecule. 
As with mode a, oscillations along this mode change the 
hydrogen bond angle and, thus, the interaction energy. 
Unlike mode a, however, the intra- molecular motions in 
this mode cause the total energy to change by a relatively 
large amount. In this case it is enough to swamp the ef- 
fect from the relatively weak energy change stemming 
from the changing hydrogen bond geometry. As a re- 
sult, all functionals find similar frequencies for this mode. 
Mode 5 corresponds to an asymmetric stretch of the hy- 
drogen acceptor molecule. This has virtually no effect 
on the hydrogen-bond geometry. The frequency of this 
mode is governed almost entirely by the intra-molecular 
strain energy in the hydrogen acceptor molecule, with 
the interaction energy playing virtually no role. Thus, 
all functionals predict very similar frequencies for this 
mode. Qualitatively, for the four modes discussed the 
spread in Etotai is proportional to the relative spread in 
frequency in Fig. [3l 

To conclude our study of the vibrational properties, 
the vdW-DF frequencies of all clusters, including the 
tetramer and the pentamer, were quantitatively com- 
pared to MP2 calculations. The results of this analysis 
are shown in Fig. [6] in the form of histograms of the er- 
rors (in percent relative to MP2), analyzing 90 modes 
in total (as in Fig. [3j modes with a frequency less than 
200 cm -1 were not included). Again, a complete list- 
ing of frequencies for these modes can be found in the 
supplementary material. In the figure a negative percent 
error simply means the calculated value was below the 
MP2 value. LDA shows a distinct bimodal distribution 
in errors. These correspond to a systematic overestima- 
tion of the frequency of low- frequency, inter-molecular 
modes and less pronounced underestimation of the high- 
frequency, intra-molecular modes relative to MP2. PBE 
has a much tighter error range but still shows the bimodal 
distribution evident in the LDA results, and for the same 
reason. As discussed above, both functionals form an ar- 
tificially strong hydrogen bond. This — together with the 
fact that inter-molecular interactions tend to enhance the 
overall restoring force for inter-molecular modes and to 
weaken it in intra-molecular modes — results in the bi- 
modal distribution in the LDA and PBE results. The 
vdW-DF plot shows an error spread that is much less se- 
vere and does not exhibit the bimodal distribution char- 
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FIG. 6. Histograms of the percent error relative to MP2 calcu- 
lations over all the frequencies of the clusters. The horizontal 
axis gives the percent error and the vertical axis gives the 
number of vibrational frequencies that had that error. 



acteristic of systematic hydrogen-bond overestimation. 
From Figs. [3] and [6] it is clear that vdW-DF system- 
atically improves upon commonly used functionals and 
demonstrates great promise to accurately reproduce vi- 
brational frequencies of water. 

IV. CONCLUSIONS 

We have used density functional theory incorporating the 
van der Waals density functional to calculate a num- 
ber of properties of small water clusters (H2 0) n with 
n = 1 — 5 and hexagonal ice Ih- Our results show ex- 
cellent agreement, both with experiment and high-level 
quantum-chemistry calculations. Since DFT is capable 
of calculating systems with a thousand atoms or more, 
this implies that accurate quantum-mechanics can now 
be applied to bulk-like water. This capability should 
rapidly advance our detailed knowledge and understand- 
ing of bulk water behavior. 

Toward this aim, the authors, working in collaboration 
with the developers of the free, open-source Quantum- 
Espresso package, 35 have implemented an efficient 48 
vdW-DF functional in the newest official release. We 
strongly encourage the testing and application of this, 
and the newer vdW-DF2 49 functional (also implemented 
in Quantum-Espresso) in water-based systems and, in- 
deed, in any system where hydrogen bonding is expected 
to play a critical role. 
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